Assignment E - Brightfield histology color deconvolution¶

  • Author: Catherine Chia and Aoming Sun
  • Teacher and TAs: Marten Postma, Aaron Lin, Aoming Sun, Catherine Chia
  • Date: 21st June, 2023

Outline of workflow¶

  1. Prerequisites:
  • Use ImageJ to crop and export images: Stain 1, Stain 2, Background, OR
  • Use ImageJ to export the RGB vectors for the same images

  1. Preprocessing

  2. Color Deconvolution

  3. Separate stains

In [19]:
#Libraries
from matplotlib import pyplot as plt, patches
import numpy as np


#Enable nice output printing features
from IPython.core.interactiveshell import InteractiveShell
InteractiveShell.ast_node_interactivity = 'last_expr_or_assign'
import warnings
warnings.filterwarnings('ignore')

#Add other libraries as you see fit
import tifffile as tif

Preprocessing steps¶

In [20]:
#Start coding here

#Import IHC image and split it to RGB
#IHC2
#img_ihc = tif.imread("IPQDA_23_ASS_E_DATA\IHC2.tif")

#H_E
img_ihc = tif.imread("IPQDA_23_ASS_E_DATA\H_E.tif")

#RGB channels
img_ihc_r = img_ihc[:, :, 0]
img_ihc_g = img_ihc[:, :, 1]
img_ihc_b = img_ihc[:, :, 2]

#Import cropped stain1, stain2 and background ROI images, OR import RGB vectors of the ROIs
#IHC2
#img_stain1 =  tif.imread("IPQDA_23_ASS_E_DATA\IHC2_stain1.tif")
#img_stain2 = tif.imread("IPQDA_23_ASS_E_DATA\IHC2_stain2.tif")
#img_background = tif.imread("IPQDA_23_ASS_E_DATA\IHC2_back.tif")
#H_E
img_stain1 =  tif.imread("IPQDA_23_ASS_E_DATA\H_E_stain1.tif")
img_stain2 = tif.imread("IPQDA_23_ASS_E_DATA\H_E_stain2.tif")
img_background = tif.imread("IPQDA_23_ASS_E_DATA\H_E_back.tif")
#End coding here
Out[20]:
array([[[251, 241, 252],
        [250, 240, 251],
        [248, 238, 247],
        ...,
        [255, 245, 255],
        [255, 247, 255],
        [255, 248, 255]],

       [[252, 242, 253],
        [247, 237, 248],
        [248, 238, 247],
        ...,
        [255, 246, 255],
        [255, 247, 255],
        [253, 243, 254]],

       [[250, 240, 251],
        [248, 238, 249],
        [249, 239, 248],
        ...,
        [255, 243, 253],
        [253, 243, 254],
        [250, 240, 251]],

       ...,

       [[249, 237, 247],
        [248, 236, 246],
        [248, 236, 246],
        ...,
        [251, 244, 251],
        [251, 244, 252],
        [252, 245, 253]],

       [[253, 241, 253],
        [252, 240, 250],
        [250, 238, 248],
        ...,
        [253, 246, 253],
        [250, 243, 251],
        [252, 245, 253]],

       [[251, 240, 254],
        [249, 239, 250],
        [251, 241, 252],
        ...,
        [255, 247, 255],
        [252, 245, 252],
        [254, 247, 254]]], dtype=uint8)
In [21]:
img_ihc.shape
Out[21]:
(2309, 3877, 3)
In [22]:
#Inspect imported IHC image
plt.title("Raw IHC image")
plt.axis('off')
plt.imshow(img_ihc)
Out[22]:
<matplotlib.image.AxesImage at 0x20280068550>

Calculate RGB mean of the images¶

In [23]:
#Start coding here

#Calculate mean of image for each RGB channels. If you use RGB vectors, assign them directly to the variables here
mean_img_stain1 = np.mean(img_stain1, axis=(0, 1))
mean_img_stain2 = np.mean(img_stain2, axis=(0, 1))
mean_img_background = np.mean(img_background, axis=(0, 1))

#End coding here

print(mean_img_stain1)
print(mean_img_stain2)
print(mean_img_background)

[249.84285714  74.72857143 117.43571429]
[66.8        13.70909091 71.57272727]
[252.23859127 243.8735119  252.17361111]

Inspect ROIs of stains and background to ensure correct stain color selection¶

In [24]:
#Convert RGB values to Hex color values for visualization
hex_img_stain1 = '#%02x%02x%02x' % tuple(mean_img_stain1.astype(int))
hex_img_stain2 = '#%02x%02x%02x' % tuple(mean_img_stain2.astype(int))
hex_img_background = '#%02x%02x%02x' % tuple(mean_img_background.astype(int))

print(hex_img_stain1)
print(hex_img_stain2)
print(hex_img_background)

#Visualization of RGB mean of cropped images
fig, axs = plt.subplots(1,3)

fig.suptitle('RGB mean of cropped images')

rectangle_stain1 = patches.Rectangle((0, 0), 1, 1, facecolor=hex_img_stain1)
rectangle_stain2 = patches.Rectangle((0, 0), 1, 1, facecolor=hex_img_stain2)
rectangle_background = patches.Rectangle((0, 0), 1, 1, facecolor=hex_img_background)

axs[0].add_patch(rectangle_stain1)
axs[1].add_patch(rectangle_stain2)
axs[2].add_patch(rectangle_background)
axs[0].set_title('Stain 1')
axs[1].set_title('Stain 2')
axs[2].set_title('Background')

axs[0].axis('off')
axs[1].axis('off')
axs[2].axis('off')
plt.show()

#f94a75
#420d47
#fcf3fc

Color Deconvolution¶

Calculate transmittance, T and convert it to absorbances, OD according to Beer–Lambert law¶

In [25]:
#Calculate transmittances, T for each stain
T_stain1 = mean_img_background/mean_img_stain1
T_stain2 = mean_img_background/mean_img_stain2
OD_stain1 = np.log10(T_stain1)
OD_stain2 = np.log10(T_stain2)

print(OD_stain1)
print(OD_stain2)

[0.00414459 0.51367795 0.33189944]
[0.57703507 1.25015598 0.54695207]

Normalize the absorbances to vector lengths¶

In [26]:
#Start coding here
#Normalize the absorbances
#Normalize by the squared of the sum of the vector values to the power of 2
OD_stain1_norm = OD_stain1 / np.max(OD_stain1)
OD_stain2_norm = OD_stain2 / np.max(OD_stain2)
#End coding here

print(OD_stain1_norm)
print(OD_stain2_norm)

[0.00806847 1.         0.64612359]
[0.46157046 1.         0.43750706]

Form a deconvolution matrix¶

In [27]:
#Start coding here

#Combine OD_stain1_norm and OD_stain2_norm to form a normalized OD matrix M
M = np.column_stack((OD_stain1_norm, OD_stain2_norm))

#Calculate the deconvolution matrix according to Linear regression
MT = np.transpose(M)

MT_M = np.dot(MT, M)

inversed_MT_M = np.linalg.inv(MT_M)

D = np.dot(inversed_MT_M, MT)
D=D.T

#End coding here

print("M")
print(M)
print("M transposed")
print(MT)
print("Inversed M transposed multiplied with M")
print(inversed_MT_M)
print("Deconvolution matrix, D")
print(D)

M
[[0.00806847 0.46157046]
 [1.         1.        ]
 [0.64612359 0.43750706]]
M transposed
[[0.00806847 1.         0.64612359]
 [0.46157046 1.         0.43750706]]
Inversed M transposed multiplied with M
[[ 4.17951775 -3.82820822]
 [-3.82820822  4.21844559]]
Deconvolution matrix, D
[[-1.73326552  1.9162221 ]
 [ 0.35130953  0.39023737]
 [ 1.02561688 -0.6278959 ]]

Calculate the coefficient for each stain¶

In [28]:
#Convert pixel intensity to transmittance to absorbance according to Beer-Lambert Law on the IHC image
#Calculate the transmittance
T_img_ihc = mean_img_background/img_ihc

#Because of the logarithmic function in the next step, we assign all transmittance value less than 1 to 1 
T_img_ihc[T_img_ihc<1] = 1
In [29]:
#Start coding here

#Calculate the absorbance
OD_img_ihc = np.log10(T_img_ihc)

#Coefficient matrix
coeffs = np.dot(OD_img_ihc, D)

#Extracting the individual coefficients from the coefficient matrix
#Which are essentially the orthogonal representation of the stains of the IHC image
coeff_stain1 = coeffs[:,:, 0]
coeff_stain2 = coeffs[:,:, 1]


#End coding here

print(coeff_stain1.shape)
print(coeff_stain2.shape)

(2309, 3877)
(2309, 3877)

Separate stains¶

Multiply the coefficients with the stain absorbance to get the image absorbance per stain¶

In [30]:
OD_img_ihc
Out[30]:
array([[[2.88995293e-02, 4.52666181e-01, 2.74594839e-01],
        [2.70631862e-02, 4.28123240e-01, 2.61820551e-01],
        [1.79961663e-02, 3.91529438e-01, 2.40331635e-01],
        ...,
        [0.00000000e+00, 0.00000000e+00, 0.00000000e+00],
        [5.61218515e-03, 0.00000000e+00, 5.50029046e-03],
        [0.00000000e+00, 0.00000000e+00, 0.00000000e+00]],

       [[5.74192586e-02, 4.57745707e-01, 2.71365869e-01],
        [3.81995524e-02, 3.91529438e-01, 2.31437922e-01],
        [1.79961663e-02, 3.37946610e-01, 2.03042551e-01],
        ...,
        [0.00000000e+00, 0.00000000e+00, 0.00000000e+00],
        [2.13781077e-03, 0.00000000e+00, 2.99096775e-04],
        [0.00000000e+00, 0.00000000e+00, 0.00000000e+00]],

       [[5.93888514e-02, 4.37774626e-01, 2.58684837e-01],
        [4.00836962e-02, 3.57780855e-01, 2.14178917e-01],
        [2.16002905e-02, 3.11617671e-01, 1.86855790e-01],
        ...,
        [0.00000000e+00, 0.00000000e+00, 0.00000000e+00],
        [0.00000000e+00, 0.00000000e+00, 0.00000000e+00],
        [0.00000000e+00, 0.00000000e+00, 0.00000000e+00]],

       ...,

       [[3.87152358e-03, 3.87164632e-01, 2.37346782e-01],
        [3.87152358e-03, 3.70131293e-01, 2.22722690e-01],
        [4.10991469e-04, 3.45771947e-01, 2.03042551e-01],
        ...,
        [5.61218515e-03, 0.00000000e+00, 0.00000000e+00],
        [1.08764251e-02, 0.00000000e+00, 3.75962888e-03],
        [4.10991469e-04, 0.00000000e+00, 0.00000000e+00]],

       [[5.61218515e-03, 3.78564461e-01, 2.34382303e-01],
        [7.35985142e-03, 3.65975333e-01, 2.22722690e-01],
        [5.61218515e-03, 3.49738134e-01, 2.08575039e-01],
        ...,
        [3.87152358e-03, 0.00000000e+00, 0.00000000e+00],
        [5.61218515e-03, 0.00000000e+00, 0.00000000e+00],
        [3.87152358e-03, 0.00000000e+00, 0.00000000e+00]],

       [[7.35985142e-03, 3.78564461e-01, 2.34382303e-01],
        [7.35985142e-03, 3.61858767e-01, 2.19856050e-01],
        [9.11457899e-03, 3.45771947e-01, 2.08575039e-01],
        ...,
        [3.87152358e-03, 0.00000000e+00, 0.00000000e+00],
        [2.13781077e-03, 0.00000000e+00, 0.00000000e+00],
        [4.10991469e-04, 0.00000000e+00, 0.00000000e+00]]])
In [31]:
#Initialize the image absorbance container per stain
OD_img_ihc_stain1 = np.zeros((img_ihc.shape[0], img_ihc.shape[1], img_ihc.shape[2]))
OD_img_ihc_stain2 = np.zeros((img_ihc.shape[0], img_ihc.shape[1], img_ihc.shape[2]))
#Start coding here
print(OD_img_ihc_stain1.shape)
print(coeff_stain1.shape)
#Multiply the coefficients with the stain absorbance per stain. Do it independently for each RGB layer"""
"""for i in range(3):
    OD_img_ihc_stain1[:, :, i] = OD_img_ihc[:, :, i] * coeff_stain1
    OD_img_ihc_stain2[:, :, i] = OD_img_ihc[:, :, i] * coeff_stain2"""
OD_img_ihc_stain1 = np.reshape(coeff_stain1, img_ihc_r.shape)[:, :, np.newaxis] * OD_stain1
OD_img_ihc_stain2 = np.reshape(coeff_stain2, img_ihc_r.shape)[:, :, np.newaxis] * OD_stain2


#End coding here

(2309, 3877, 3)
(2309, 3877)
Out[31]:
array([[[0.03439603, 0.07451955, 0.03260283],
        [0.03146729, 0.0681744 , 0.02982678],
        [0.02098724, 0.0454692 , 0.0198931 ],
        ...,
        [0.        , 0.        , 0.        ],
        [0.00421269, 0.00912687, 0.00399307],
        [0.        , 0.        , 0.        ]],

       [[0.06824479, 0.14785347, 0.06468694],
        [0.04654903, 0.10084924, 0.04412225],
        [0.02243193, 0.04859914, 0.02126247],
        ...,
        [0.        , 0.        , 0.        ],
        [0.00225547, 0.00488651, 0.00213788],
        [0.        , 0.        , 0.        ]],

       [[0.07052009, 0.15278294, 0.06684362],
        [0.0472861 , 0.10244611, 0.0448209 ],
        [0.02635309, 0.05709441, 0.02497921],
        ...,
        [0.        , 0.        , 0.        ],
        [0.        , 0.        , 0.        ],
        [0.        , 0.        , 0.        ]],

       ...,

       [[0.00546783, 0.01184615, 0.00518277],
        [0.00693083, 0.01501575, 0.0065695 ],
        [0.00474964, 0.01029016, 0.00450202],
        ...,
        [0.00620555, 0.01344442, 0.00588203],
        [0.01066418, 0.02310412, 0.01010822],
        [0.00045444, 0.00098456, 0.00043075]],

       [[0.00653002, 0.0141474 , 0.00618959],
        [0.00985212, 0.02134479, 0.0093385 ],
        [0.00938933, 0.02034213, 0.00889983],
        ...,
        [0.00428085, 0.00927453, 0.00405767],
        [0.00620555, 0.01344442, 0.00588203],
        [0.00428085, 0.00927453, 0.00405767]],

       [[0.00846246, 0.01833407, 0.00802128],
        [0.00996379, 0.02158671, 0.00944434],
        [0.01236891, 0.02679745, 0.01172407],
        ...,
        [0.00428085, 0.00927453, 0.00405767],
        [0.00236384, 0.00512129, 0.0022406 ],
        [0.00045444, 0.00098456, 0.00043075]]])

Convert the image absorbance to image transmittance¶

In [32]:
#Convert absorbance to transmittance
T_img_ihc_stain1 = 10**(-OD_img_ihc_stain1)
T_img_ihc_stain2 = 10**(-OD_img_ihc_stain2)
Out[32]:
array([[[0.92385534, 0.84232646, 0.92767781],
        [0.93010657, 0.85472341, 0.9336266 ],
        [0.95282416, 0.90059763, 0.95522769],
        ...,
        [1.        , 1.        , 1.        ],
        [0.99034681, 0.97920389, 0.99084776],
        [1.        , 1.        , 1.        ]],

       [[0.85458488, 0.71145351, 0.86161462],
        [0.89836117, 0.79277649, 0.90339513],
        [0.94965984, 0.89413039, 0.95222051],
        ...,
        [1.        , 1.        , 1.        ],
        [0.99482006, 0.98881146, 0.99508944],
        [1.        , 1.        , 1.        ]],

       [[0.85011936, 0.7034238 , 0.85734651],
        [0.89683779, 0.78986685, 0.90194302],
        [0.94112412, 0.87681019, 0.94410607],
        ...,
        [1.        , 1.        , 1.        ],
        [1.        , 1.        , 1.        ],
        [1.        , 1.        , 1.        ]],

       ...,

       [[0.98748877, 0.97309188, 0.98813715],
        [0.98416785, 0.96601585, 0.98498701],
        [0.98912315, 0.97658452, 0.98968726],
        ...,
        [0.9858128 , 0.96951734, 0.98654743],
        [0.97574384, 0.94819111, 0.97699375],
        [0.99895415, 0.99773553, 0.99900865]],

       [[0.98507654, 0.96794928, 0.98584903],
        [0.97757003, 0.95204004, 0.97872685],
        [0.97861231, 0.95424056, 0.97971594],
        ...,
        [0.9901914 , 0.97887102, 0.99070037],
        [0.9858128 , 0.96951734, 0.98654743],
        [0.9901914 , 0.97887102, 0.99070037]],

       [[0.98070307, 0.95866292, 0.98169983],
        [0.97731871, 0.95150985, 0.97848836],
        [0.97192128, 0.94016169, 0.97336545],
        ...,
        [0.9901914 , 0.97887102, 0.99070037],
        [0.99457185, 0.98827705, 0.99485411],
        [0.99895415, 0.99773553, 0.99900865]]])

Clip each layer in the image transmittance to values between 0 and 1, preparing for conversion to values between 0 and 255 later¶

In [33]:
#Clip each layer to 0,1
T_img_ihc_stain1[T_img_ihc_stain1>1] = 1
T_img_ihc_stain2[T_img_ihc_stain2>1] = 1
T_img_ihc_stain1[T_img_ihc_stain1<0] = 0
T_img_ihc_stain2[T_img_ihc_stain2<0] = 0

Convert the image transmittance to values between 0 and 255 (integers), so that plotting is possible¶

In [34]:
#Start coding here

T_img_ihc_stain1_norm = (T_img_ihc_stain1 * 255).astype(np.uint8)

T_img_ihc_stain2_norm = (T_img_ihc_stain2 * 255).astype(np.uint8)
#End coding here
Out[34]:
array([[[235, 214, 236],
        [237, 217, 238],
        [242, 229, 243],
        ...,
        [255, 255, 255],
        [252, 249, 252],
        [255, 255, 255]],

       [[217, 181, 219],
        [229, 202, 230],
        [242, 228, 242],
        ...,
        [255, 255, 255],
        [253, 252, 253],
        [255, 255, 255]],

       [[216, 179, 218],
        [228, 201, 229],
        [239, 223, 240],
        ...,
        [255, 255, 255],
        [255, 255, 255],
        [255, 255, 255]],

       ...,

       [[251, 248, 251],
        [250, 246, 251],
        [252, 249, 252],
        ...,
        [251, 247, 251],
        [248, 241, 249],
        [254, 254, 254]],

       [[251, 246, 251],
        [249, 242, 249],
        [249, 243, 249],
        ...,
        [252, 249, 252],
        [251, 247, 251],
        [252, 249, 252]],

       [[250, 244, 250],
        [249, 242, 249],
        [247, 239, 248],
        ...,
        [252, 249, 252],
        [253, 252, 253],
        [254, 254, 254]]], dtype=uint8)

Visualize deconvolved images¶

In [35]:
#Display deconvolved image for stain 1
fig = plt.figure(dpi=300)
plt.title("Deconvolved image for stain 1")
plt.axis('off')
plt.imshow(T_img_ihc_stain1_norm)
fig.savefig('T_img_ihc_stain1_norm.tif')
In [36]:
#Display and export deconvolved image for stain 1
fig = plt.figure(dpi=300)
plt.title("Deconvolved image for stain 2")
plt.axis('off')
plt.imshow(T_img_ihc_stain2_norm)
fig.savefig('T_img_ihc_stain2_norm.tif')
In [36]: